arXiv:1507.01751vl [cond-mat.mtrl-sci] 7 Jul 2015 


Ab initio calculation of the shock Hugoniot of bulk silicon 


Oliver StricksoiQ 

Laboratory for Scientific Computing, Cavendish Laboratory, University of Cambridge, 

J. J. Thomson Avenue, Cambridge CBS OHE, United Kingdom 

Emilio Artacho 

Theory of Condensed Matter, Cavendish Laboratory, University of Cambridge, 

J. J. Thomson Avenue, Cambridge CBS OHE, United Kingdom 
CIC Nanogune and DIPC, Tolosa Hiribidea 76, 20018 San Sebastian, Spain and 
Basque Foundation for Science Ikerbasque, Bilbao, Spain 
(Dated: July 2015) 

We describe a simple annealing procedure to obtain the Hugoniot locus (states accessible by a shock wave) for 
a given material in a computationally efficient manner. We apply this method to determine the Hugoniot locus 
in bulk silicon from ab initio molecular dynamics with forces from density-functional theory, up to 70 GPa. The 
fact that shock waves can split into multiple waves due to phase transitions or yielding is taken into account here 
hy specifying the strength of any preceding waves explicitly based on their yield shain. Points corresponding 
to uniaxial elastic compression along three crystal axes and a number of post-shock phases are given, includ¬ 
ing a plastically-yielded state, approximated hy an isotropic stress configuration following an elastic wave of 
predetermined strength. The results compare well to existing experimental data for shocked silicon. 


I. INTRODUCTION 

Shock waves are used extensively to study matter at con¬ 
ditions of extreme pressure and temperature, and have been 
used to obtain some of the highest laboratory-attained pres¬ 
sures. They are useful for equation of state determination and 
are important dynamic phenomena in their own right, arising 
in aerodynamics,!^ reactive flovJ^land high-speed impact.® 

Simulations of shock waves have a long history.l^ Direct 
simulations using empirical potentials are now feasible on a 
multi-billion atom scale on present hardware, which is large 
enough to observe detailed mechanisms of yield, plastic flow 
and shock interaction with nanostructures, directly.® Work 
with empirical potentials can give important insight and un¬ 
derstanding, but a need for first-principles methods such as 
Density Functional Theory (dft) exists in providing predic¬ 
tive power and accuracy. These methods must use more mod¬ 
est system sizes, of hundreds or thousands of atoms in the case 
of DFT. 

Silicon has a rich phase diagram, with metallic dense 
phases rather different in character to the ambient diamond 
phase, making it an interesting and challenging object of sim¬ 
ulation. In total, eleven stable or metastable phases of silicon 
are currently known.® Shock experiments have provided im¬ 
portant data for constructing the phase diagram. The phase 
transition in silicon from the cubic diamond structure to the 
beta-tin structure, occurring at 12 GPa at room temperature, 
and undergoing a reduction in volume of 20%, has been 
well est ablis hed by static loading experiments from the 1960s 
onward.lSnSI Evidence of at least one phase transition at sim¬ 
ilar pressures was then observed in shock-wave experiments, 
starting with Pavlovskii .E! 

If a shock wave is strong enough to cause a material to yield 
plastically or undergo a phase transition, the wave can split 
into two or more separate shock waves, and this has long been 
observed and understood.® In this situation, the last shock 
takes the material to its final state, but the preceding shocks 


take the material to a cusp on the pressure-volume Hugoniot 
locus caused by a transition: either the Hugoniot elastic limit 
or the onset pressure of a phase transition. In silicon. Gust and 
RoyceEnH found a three-wave structure for samples shocked 
in the (100) crystal direction and a four-wave structure when 
shocked in the (110) or (111) directions. In the latter cases, 
these waves were attributed to: an initial elastic precursor to 
the Hugoniot elastic limit of 5.5 GPa, followed by waves cor¬ 
responding to a state of plastic yield and two successive phase 
transitions at 10 GPa and 13 GPa. Along (100), the higher 
elastic limit of 9 GPa obscures the first transition wave, and a 
single wave takes the material simultaneously to a new phase 
and to a state of hydrostatic stress. 

The work of Goto et a/. !SI largely confirmed the findings of 
Gust and Royce ,121 although they observed a three-wave struc¬ 
ture, regardless of crystal orientation, consistent with only a 
single phase transition at 13 GPa. Above the Hugoniot elastic 
limit, shock compression was found to result in a hydrostatic 
stress configuration, due to the complete loss of strength in the 
material. 

More recently, and con trary to the earlier experimental 
work, Turneaure and Gupta^i^^ reported a single phase tran¬ 
sition that is complete by 15.9 GPa. Shocks to these pressures 
show a much greater volume compression than the points at¬ 
tributed to an extended mixed-phase region by both Gust and 
Royce^i^ and Goto et Here the phase transition is not 
complete until at least 30 GPa. This discrepancy is explained 
by Turneaure and Gupta^i^ as arising from the reflection of the 
first two shock waves propagating back into the material be¬ 
fore the arrival of the third wave, and altering the peak state of 
the earlier experiments. They avoid this eventuality by back¬ 
ing the silicon with a window made from lithium fluoride, a 
material with a good impedance match to silicon. 

The Imma phase of silicon is found intermediate between 
the beta-tin and simple hexagonal phases, and is stable be¬ 
tween 13 GPa and 15 GPa at room temperature .^3 Theoreti¬ 
cally, the energy and volume of these three phases are close.® 
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A recent simulation of directly shocked silicon using an em¬ 
pirical potential^ found a phase transition to an Imma phase 
with a modification of the Tersoff potentiafSEHl of Erhart and 
Albe P 

In this paper, we give the Hugoniot loci according to Den¬ 
sity Functional Theory for several pure phases of silicon, 
including cubic diamond under elastic compression along 
(100), (110) and (111), a hydrostat (resulting from either 
a single shock or a split-shock structure), beta-tin, simple 
hexagonal and the liquid, and report shock temperatures for 
these states. 

Several approaches can been taken for the determination 
of a Hugoniot locus from molecular dynamics. The most 
straightforward, but computationally the most demanding, is 
to simulate a slab of atoms struck by an impactor directly, 
measuring the speed of any shock waves and post-shock aver¬ 
age particle velocities as they arise from the simulation. From 
the Hugoniot relations, these velocities can be converted to a 
relationship between pressure and volume compression. For 
empirical potentials, a local stress is conveniently available, 
so this could also be taken directly from the simulation. This 
is the approach taken by, e.g. Kadau et 

It is simple to check that a given equilibrium state lies on or 
close to the (single-shock) Hugoniot locus, which amounts to 
satisfying the Hugoniot relation 

£-£o = -^(a^^ + CTo”)(yo-V), (1) 

where E is the internal energy, V is the specific volume and 
(7^^ is the stress in the direction of the shock (and can be re¬ 
placed with the pressure p in a hydrostatic situation). The 
zero-subscripted variables are for the pre-shocked state. Other 
(equivalent) Hugoniot relations exist between any three of; in¬ 
ternal energy, pressure, volume, shock velocity and particle 
velocity. It is therefore sufficient to sample several points that 
are chosen to bracket the Hugoniot locus, and the Hugoniot 
state then approximated by interpolation, or solved for itera¬ 
tively. The former is the approach taken by Bonev et 
shocked deuterium. 

Swift et aZ.ES constructed a polymorphic equation of state 
for silicon, incorporating DFT simulations of the cubic dia¬ 
mond and j3-Sn phases, with the lattice-thermal contribution 
approximated by quasiharmonic phonons. The equation of 
state was constructed with a particular focus on simulating 
shock waves. The full equation of state was sampled and 
the Hugoniot locus could therefore be extracted as a one¬ 
dimensional path through it. The phase boundary and mixed 
phase region along the Hugoniot were found explicitly by 
minimizing the Helmholtz free energy computed from the 
quasiharmonic phonon approximation. 

Alternatively, a Hugoniot state can be determined dynam¬ 
ically from within a single molecular dynamics simulation 
by some modified dynamics to constrain the state to sat¬ 
isfy eq. fTb. This is the approach taken by the Hugoniostat 
methodand the technique of Reed et The former 
simulations use modified Nose-Hoover dynamics while the 
latter uses coupled dynamics of the atoms and simulation cell, 
whose Lagrangian involves the computed instantaneous shock 


TABLE I. Basis parameters for silicon, according to the soft- 
confinement scheme of Junquera et fll.EI. For the purposes of basis 
generation, an effective ionic charge of -0.46 was used, which was 
also variationally optimized. The cutoff radii of the first and second 
zeta functions are r(^i) and r{^2), and r, is the confinement poten¬ 
tial’s internal radius. Vq is the soft-confinement prefactor. 


n 

/ 

G (up) 

r(^l)(flp) 

r(C 2 ) (ao) 

Vo (Ry) 

3 

0 

4.97 

7.00 

4.38 

15.43 

3 

1 

3.83 

7.00 

4.09 

4.70 

3 

2 

0.03 

4.55 

- 

11.97 


speed, and varies the simulation cell uniaxially. One aim of 
these dynamics is to work on timescales comparable to shock- 
passage times, without the overhead of dealing with a direct 
non-equilibrium simulation. 

If we are interested only in the final post-shock state, and 
are not interested in the (modified) dynamics while the con¬ 
straint is being applied, we are free to use a method based 
on simple velocity rescaling, analogous to the procedure of 
Berendsen,l2Hl which is what we propose here due to its in¬ 
creased efficiency in reaching the final state. 


II. COMPUTATIONAL METHOD 
A. Density Functional Theory 

The ah initio MD simulations described here were per¬ 
formed with the Siesta method and implementation of Den¬ 
sity Functional Theory,!^ using the Perdew et GGA func¬ 
tional. 

The core electrons were described with a Troullier-Martins 
norm-conserving pseudopotentiai^with a matching radius in 
each angular momentum channel of 1.89 ap. The valence 
electrons were described with a basis of numerical atomic or¬ 
bitals of double-^ polarized typd^H (representing 13 orbitals 
per atom). The basis was generated by fixing the longest or¬ 
bital cutoffs at 7.0 ap and variationally optimising the other 
parameters in bulk diamond-phase silicon—the final basis pa¬ 
rameters are given in table 

The mesh used for integrals in real-space was well con¬ 
verged at a grid cutoff of 100 Ry. The dense phases of silicon 
required several k-points to converge in energy, and in partic¬ 
ular, for the cold compression curves of the various phases to 
converge in energy relative to one another. A 4^ Monkhorst- 
Pack grid of points was used on the 64 atom simulations, to 
give an effective cutoff length of 21 A. 

The electronic temperature used in the DFT calculations 
should be consistent with the final temperature attained after 
the annealing process described below. The consistent forces 
for the ah initio molecular dynamics are the nuclear-position 
derivatives of the electronic free-energy as defined in Mer- 
min’s DFT.ISI All of the simulations reported below are for 
an electronic temperature of 300 K, except for the two points 
with highest temperatures, for which the electronic tempera¬ 
ture was adjusted to coincide with the final (nuclear) temper- 
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ature. The effect of the electronic temperature on the reported 
quantities was found to be quite small: the maximum differ¬ 
ence in pressure for the hottest simulation between using a 
consistent electronic temperature and the initial 300 K is be¬ 
low 5%. 

The integration of the dynamics used the Born- 
Oppenheimer approximation with a timestep of 1 fs. 


B. Annealing to the Hugoniot Locus 


We use a simple annealing procedure to find the state on the 
Hugoniot corresponding to a specified longitudinal strain. A 
Berendsen thermostaP^ is used with a variable target temper¬ 
ature computed from the instantaneous difference in energy 
between the total energy of the system, and the total energy 
that would be required to satisfy the energy Hugoniot relation, 
eq. exactly, given the current instantaneous longitudinal 
stress. 

The procedure is given explicitly below. This may be com¬ 
bined with a further anneal to relax the pressure to a hydro¬ 
static configuration if desired. Optionally, the box vectors 
may be gradually ramped between two states, which is most 
useful when the starting state of the simulation and the initial 
state of the Hugoniot locus are the same. 


procedure HugoniotAnneal(£'q°‘, Vq, cto) 
compute g,Fn from atomic positions jc„ 
for all atoms do [> velocity Verlet 

Vfi -)- — (^Ffj—i -f Ffi) 


y(unc) 

'■n+l 


2m 

^ jc„ -I- df v„ 


end for 

pL 

pkin , y 


^Fnjm 


,mv„ 


atoms 2 


mv, 


)v„ 

• V„ 


> compute the stress 


> compute the target energy 
£hug ^ £tot _ 1 ((^33 (j 33) _ y ) 


^target ^ ~ 




1 + 


dt 


Trelax \ E^' 


ykin 

target ^ 


for all atoms do 


> scale the velocities 
[> correct positions based on the scaled velocities 

,(unc) , . ^.(sca) 


end for 

t i — t -f dt, n i — n-\- 1 

end procedure 


The meaning of the variables used is as follows. E denotes 
an energy (refer to the sub and superscripts), V is the unit cell 
volume, (7 is the stress, which is the sum of a kinetic term and 
the strain derivative of the total electronic energy g',x„,Vn and 
F„ are the atomic positions, velocities and forces at the nth 
timestep (‘unc’ stands for ‘uncorrected’ and ‘sea’ for scaled), 
m is the mass of a given atom, Treiax is the relaxation time, t 
and df are the current time and timestep, and anything with 
a subscript ‘0’ refers to its (time averaged) value in the un¬ 
shocked state (which may be different from the starting state 
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FIG. 1. Response of (a) internal energy and the difference with the 
Hugoniot energy computed from eq. (b) temperature and target 
temperature, and (c) pressure, to the Hugoniot anneal described here 
with a relaxation time of 100 fs. After 2000 fs, the anneal is switched 
off and the dynamics continued with Verlet integration. The response 
is averaged over 10000 independent 216 atom Stillinger-Weber sili¬ 
con systems, starting from a 2000 K liquid and annealed to the Hugo¬ 
niot locus with an initial state of 300 K and zero stress. For compar¬ 
ison, the light grey lines are taken from a single trajectory—in the 
energy plot, this is indistinguishable from the mean. 


of the simulation). 

Even though Berendsen thermo- and barostats do not repro¬ 
duce canonical statistics,ISlit is well known that they are much 
more efficient at annealing to a given equilibrium state at a de¬ 
sired temperature or pressure compared with modified dynam¬ 
ics, such as Nose-Hoover. The same applies here, compared 
to the related Hugoniostat for shocks, and this justifies their 
use, since we are interested only in the outcome of the anneal, 
not the intermediate dynamics. After the time-averaged state 
of the system closely satisfies the Hugoniot relation, the simu¬ 
lation can be restarted with Verlet dynamics to check if eq. ([T]) 
is indeed satisfied. 

Figure[T]shows the convergence in total energy, temperature 
and pressure of liquid Stillinger-Weber silicon to a state on the 
particular Hugoniot locus from an initial state of 300 K and 
zero pressure. This is an averaged result of 10000 indepen¬ 
dent simulations, each starting in the liquid phase at 2000 K. 
The relaxation time used was 100 fs. After 2000 timesteps of 
1 fs, the anneal is switched off and Verlet integration used for 
the remaining time. Note the slight relaxation of temperature 
and pressure away from their final values under the thermo- 
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FIG. 2. Pressure-volume Hugoniot loci for silicon. The solid red 
lines in the figure are the DFT results from this work (with contained 
points indicating individual simulations), with an initial pre-shocked 
state of zero pressure and 300 K, with the final state in the indicated 
phase (‘sh’ for simple hexagonal). Estimated error is less than 5% 
for the liquid and beta-tin phases, and is substantially smaller for the 
diamond phase. The symbols are experimental results from the litera¬ 
ture: o Gust and Rovcel^. o Goto et a/.G^, a Tumeaure and Gupta^^, 
A Tumeaure and Gupta^^^. The dashed lines are approximations to 
the mixed-phase portion of the Hugoniot, for cubic diamond to: liq¬ 
uid (upper blue line) and beta tin (lower green line). 


FIG. 4. Particle velocity-shock velocity Hugoniot loci for silicon. 
The DFT results (solid red lines and points) each correspond to an 
initial state of zero pressure and 300 K, with the final state in the 
indicated phase. The dashed line is for a single-shock process whose 
final state has a hydrostatic stress configuration. The meaning of 
the symbols is the same as in fig. with the blue diamonds on the 
axis the elastic and bulk wave speeds from Goto et The dotted 
base line indicates equal shock and particle velocity, below which no 
viable shock should be recorded. 




Relative volume 


FIG. 3. Pressure-volume Hugoniot loci for silicon. This is a similar 
plot to fig. 1^ with the meaning of the symbols and lines the same, 
emphasizing the small strain region of the Hugoniot locus and with 
the results of Gust and Roycel^^ omitted due to their larger variance. 


FIG. 5. Post-shock temperature as a function of volume for several 
final states. The DFT results (solid red lines) each correspond to an 
initial state of zero pressure and 300 K, with the final state in the 
indicated phase. The ‘plastic’ curve does not include the temperature 
rise due to dissipative heating. The meaning of the symbols is the 
same as in fig.|^ 


stat. For this case, it amounts to a temperature difference of 
within 1%. 


III. RESULTS 

The calculated pressure-volume and shock-velocity- 
particle-velocity Hugoniot loci for the pure phases are com¬ 


pared to results from several experiments in figs.toThe 
specific volume at zero pressure and 300 K for the PBE func¬ 
tional is 0.421 cm^/g, which is smaller than the experimental 
value of 0.431 cm^/g. The reduced volume is plotted in the 
figures; if the specific volume were plotted instead, the DFT 
results would be offset by an amount corresponding to the dif¬ 
ference in zero-pressure volume. The particle velocity-shock 
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TABLE II. Coefficients of a linear fit the shock velocity for the elastic 
waves, Us = cq +sup, for this work and two sets of experimentally- 
determined values. 

(100) (110) (111) Bulk 

CQ S Co S Cq S Co s 


(km/s) (-) 


This work 

8.38 

0.42 

9.21 

0.57 

9.34 

0.57 

6.51 1.18 

Ref.Ql 

8.42 

0.32 

9.24 

1.01 

9.39 

0.98 

- 

Ref. [35] 

8.43 

- 

9.13 

- 

9.3f] 

- 

6.4^ - 


“ Calculated from the given elastic constants and density. 


velocity Hugoniot is not greatly affected. 

The curves for the elastic shocks are computed from a 
uniaxial box deformation along the indicated direction. The 
‘plastic’ curve is for a split shock, with an elastic precursor 
to 6 GPa, taking the material to a hydrostatic stress configura¬ 
tion; this supposes that the material has no residual strength. 
The hydrostat in fig. is for an unphysical shock process 
that relaxes the material to hydrostatic stress behind a single, 
unsplit shock wave. This permits comparison with the bulk 
speed of sound (the shock velocity for this wave should ex¬ 
trapolate to the bulk speed of sound at zero particle velocity.) 

When comparing the hydrostat and the ‘plastic’ curve to 
the yielded phase, we assume that the yielding serves only to 
remove the deviatoric stress, and that the bulk response of the 
material is unaffected. We neglect the dissipative heating due 
to this effect. 

The agreement with the experimental data for the elastic 
and plastic shocks is good, with the compressibility along 
(100), (110) matching well in value and (111) showing the 
correct trend (although underestimating the value). The close 
match between the experimental plastic shock pressures and 
the hydrostatic plastic shock calculated here supports the ob¬ 
servation that the material loses all of its strength after yield. 

The particle and shock velocities in fig. are computed 
from the computed pressure and volume points using the 
Hugoniot relations 

u], = {p-po){vo-v) (2) 

Us =vl{p-po)l{vo-v), (3) 

where vq and v are the initial and final specific volumes. A 
linear fit to the elastic part of the shock-velocity-particle- 
velocity Hugoniot has coefficients given in table The ex¬ 
trapolated value of the bulk sound speed of 6.51 kms^^ agrees 
very well with the value of 6 .48 k ms~^ calculated from the 
second order elastic constants.G^EU 

The pSn and simple hexagonal curves each correspond 
to a three wave split shock structure, behind an elastic wave 
to the experimental elastic limit of 6 GPa and a secondary 
wave to the experimental location of the phase transition at 
13.8 GPa. For both of these waves, the computed volume 
for the (100) direction was used for the post-shock state. In 
general, it is quite insensitive to the precise location of the 
wave split, particularly for the elastic case, since the contri¬ 
bution to the energy change is much smaller than the 20% 
volume reduction across the phase change. The final stress 


was hydrostatic. Since the c/a-ratio is free in the )3-Sn and 
simple hexagonal structures, an additional relaxation step was 
used on the simulation box to impose a hydrostatic distribu¬ 
tion of stress while simultaneously annealing to the Hugoniot. 
The j3-Sn and simple hexagonal curves are close in pressure, 
temperature and shock velocity, with the experimental values 
closest to the simple hexagonal DFT Hugoniot. The computed 
pressures and temperatures of these points put them in stable 
region for the simple hexagonal structure on the silicon phase 
diagram.l^SI 

Part of the liquid Hugoniot corresponds to a three-wave 
shock structure, with the third wave reaching the final liquid 
state, behind a secondary wave to the onset of the melting 
transition and an elastic precursor wave. For the highest pres¬ 
sures, where the final wave has a velocity greater than that 
of the secondary wave of 6.83 km/s, it instead corresponds 
to a two-wave structure (behind only the elastic precursor). 
The largest shock pressures closely match the calculated liq¬ 
uid Hugoniot, with the simulated liquid being systematically 
slightly too stiff. 

The predicted post-shock temperatures (given in fig. in¬ 
dicate that these highest pressure points are likely to be liquid 
phase. The sixfold coordinated liquid lies close in p-v to the 
Hugoniot for the beta-tin phase, and so this phase transition 
does not exhibit the large mixed phase region as for the dia¬ 
mond to dense-phase silicon. 

A. The Phase Transition 

There is a considerable range of relative volume between 
the Hugoniot loci of the pure phases shown in fig. The ex¬ 
perimentally measured points in this region have a final state 
that is a mixture of two phases. Points on the mixed-phase 
region of the Hugoniot are on the intersection of the phase 
boundary for the two phases, as well as satisfying eq. Q. 

Similar to the plastic shock, a pressure-volume Hugoniot 
is convex at the onset of a mixed phase region; if the change 
in slope is great enough, this causes the shock to split into a 
wave taking the material to the pressure at the onset of the 
phase transition, and a slower wave taking the material to its 
final state, which is a coexistence of the two phases. 

The Hugoniot locus through the mixed phase region can 
be constructed by considering the jump condition in enthalpy 
across the shock from the point (‘1’) at the onset of the transi¬ 
tion to a point (‘2’) on the mixed Hugoniot 

H 2 —Hi = El —Ei+P 2 V 2 —piVi, (4) 

and on substituting eq. 0 for the jump in internal energy, this 
reduces to 

H2-Hi = ]^{p2-pi)iy2+Vi). (5) 

The latent heat L of the phase transition results in a change in 
enthalpy, written according to the Clausius-Clapeyron equa¬ 
tion as 

XL=-T%{Vi-V2), 


(6) 
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TABLE III. Summary of values used at the onset of the cubic dia¬ 
mond to liquid phase transition. The phase line is as obtained by 
the experiment of Kubo et al^^. The other values are from Hull®, 
with a and Cp at 1600 K and ambient pressure, and j6 at 298 K and 
13.8 GPa._ 

T (K) dT/dP (KGPa-') a (K^') /3 (GPa-') Cp (Jg- * ') 

1683 6T4 4.5 x lO^'’ 0.024 TO 


where A is the mass fraction of the second phase and the 
derivative is along the phase line. 

Since the mixed region is not at constant pressure, there is 
an additional contribution to the enthalpy change from the dif¬ 
ference in pressure and volume between the onset of the tran¬ 
sition and the post-shock state. This leads to a linearized equa¬ 
tion relating the pressure and volume changes on the phase- 
transition shock,® 


P2-P\ = (Vi -y2)x 

-1 


(7) 

where p is the isothermal compressibility, a is the volumetric 
thermal expansion coefficient and Cp is the heat capacity at 
constant pressure. The derivative dT/dp is once again along 
the phase boundary. 

We require knowledge of the onset of the transition in the 
p-V plane, which is not available from the single phase sim¬ 
ulations alone (the simulated materials are capable of be¬ 
ing substantially superheated or supercooled). This could be 
obtained from the point where the Hugoniot cuts the phase 
boundary obtained by some other method. 

We consider here two possible phase transitions starting 
from silicon in the cubic diamond structure; to a liquid, and 
to the beta-tin structure. In addition, we assume that the onset 
of either transition occurs at 13.8 GPa, close to the observed 
experimental value. Th^hase lines are experimental values, 
obtained by Kubo et a/.® This gives the two dashed lines ap¬ 
pearing in fig. The lower, green dashed line for diamond 
structure to beta-tin is nearly at constant pressure, since its 
slope is dominated by the steep phase-line of the transition,® 
dT/dp = —1426 K/GPa. This is consistent with the experi¬ 
ment of Tumeaure and Gupta .® The upper, blue dashed line 
for melting the diamond structure is influenced most strongly 
by the compressibility p of the cubic diamond phase at the 


n , 1 ^ \ dr C„ /dT 

W+l-(V,-V,)-2aV,j- + ^(^- 


pressure and temperature of the onset. Representative litera¬ 
ture values for the constants appearing in the above expression 
for the liquid are summarized in table This line underes¬ 
timates the experimentally observed slope seen by Gust and 
Royce® and Goto et a/.® While the simulated temperature 
at this pressure is much too low for melting, the simulations of 
the ‘plastically-yielded’ state do not include dissipative heat¬ 
ing and this could cause a considerable temperature rise above 
those reported in fig. 


IV. CONCLUSION 

In conclusion, we have described a simple annealing 
method and shown that it may be used to obtain a state on 
the Hugoniot locus of a pure phase of a material with sev¬ 
eral condensed phases efficiently, from first-principles. An 
approximation relying on the slope of the phase boundary can 
be used to obtain the part of the Hugoniot corresponding to 
coexistence between two phases. 

In the case of silicon, the results computed using this proce¬ 
dure with the forces described using density functional theory 
match existing experimental data very well for pressures up 
to 60 GPa, the limit of available experimental data. We have 
provided a prediction of the shock temperatures of silicon over 
this pressure range. This study supports the conclusions of the 
experimental work in general, that silicon after yield supports 
no deviatoric stress, and of Turneaure and Guptathat the 
first observed phase transition along the shock locus is likely 
to be to simple hexagonal. 


ACKNOWLEDGMENTS 


This research was supported with funding from Orica 
Ltd. and the following grants: MINECO-Spain’s Plan Na- 
cional Grant No. FIS2012-37549-C05-01, Basque Govern¬ 
ment Grant PI2014-105 CIC07 2014-2016, EU Grant “Elec- 
tronStopping” in the Marie Curie CIG Program. Part of 
this work was performed using the Darwin Supercomputer of 
the University of Cambridge High Performance Computing 
Service (http: //www.hpc. carni. ac.uk/), provided by Dell 
Inc. using Strategic Research Infrastructure Funding from the 
Higher Education Funding Council for England and fund¬ 
ing from the Science and Technology Facilities Council. We 
thank Alan Minchinton, Richard Needs, Nikos Nikiforakis, 
Stephen Walley and David Williamson for useful input and 
discussions. 


* ots22@cam.ac.uk 

' D. S. Dolling, AIAA journal 39, 1517 (2001) 

2 D. D. Dlott, Ann. Rev. Phys. Chem. 62, 575 (2011) 

3 G. E. Duvall and R. A. Graham, Rev. Mod. Phys. 49, 523 (1977) 
J. Asay and M. Shahinpoor, eds., High-Pressure Shock Com¬ 
pression of Solids, Shock Wave and High Pressure Phenomena 
(Springer-Verlag New York, 1993). 


5 B. Holian, Shock Waves 13, 489 (2004) 

® K. Kadau, T. C. Germann, and P. S. Lomdahl, Int J Mod Phys C 
17, 1755 (2006) 

' A. Shekhar, K.-i. Nomura, R. K. Kalia, A. Nakano, and 
P. Vashishta, Phys. Rev. Lett. Ill, 184503 (2013) 

* A. Mujica, A. Rubio, A. Munoz, and R. J. Needs, Rev. Mod. 
Phys. 75, 863 (2003) 












7 


^ S. Minomura and H. Drickamer, J. Phys. Chem. Solids 23, 451 
(1962) 

J. C. Jamieson, Science 139, 762 (1963) 

" M. Pavlovskii, Sov. Phys. Solid State 9, 2514 (1968). 

W. Gust and E. Royce, Dynamic yield strengths of light armor ma¬ 
terials., Tech. Rep. (California Univ., Livermore. Lawrence Radi¬ 
ation Lab., 1970). 

W. Gust and E. Royce, J. Appl. Phys. 42, 1897 (1971) 

T. Goto, T. Sato, and Y. Syono, Jpn. J. Appl. Phys. 21, L369 
(1982) 

S. J. Turneaure and Y. Gupta, Appl. Phys. Lett. 91,201913 (2007) 
S. J. Turneaure and Y. Gupta, Appl. Phys. Lett. 90,051905 (2007) 
M. 1. McMahon, R. J. Nelmes, N. G. Wright, and D. R. Allan, 
Phys. Rev. B 50, 739 (1994) 

G. Mogni, A. Higginbotham, K. Gaal-Nagy, N. Park, and J. S. 
Wark, Phys. Rev. B 89, 064104 (2014) 

J. Tersoff, Phys. Rev. Lett. 56, 632 (1986) 

20 J. Tersoff, Phys. Rev. B 38, 9902 (1988) 

21 P. Erhart and K. Albe, Phys. Rev. B 71, 035211 (2005) 

22 K. Kadau, T. C. Germann, P. S. Lomdahl, R. C. Albers, J. S. Wark, 
A. Higginbotham, and B. L. Holian, Phys. Rev. Lett. 98, 135701 
(2007) 

2^ S. A. Bonev, B. Militzer, and G. Galli, Phys. Rev. B 69, 014101 
(2004) 

2^ D. C. Swift, G. J. Ackland, A. Hauer, and G. A. Kyrala, Phys. 
Rev. B 64, 214107 (2001) 


2^ J.-B. Maillet, M. Mareschal, L. Soulard, R. Ravelo, P. S. Lom¬ 
dahl, T. C. Germann, and B. L. Holian, Phys. Rev. E 63, 016121 
( 2000 ) 

20 R. Ravelo, B. L. Holian, T. C. Germann, and P. S. Lomdahl, Phys. 
Rev. B 70, 014103 (2004) 

2' E. J. Reed, L. E. Fried, and J. D. Joannopoulos, Phys. Rev. Lett. 
90, 235503 (2003) 

2*1 H. J. C. Berendsen, J. P. M. Postma, W. E. van Gunsteren, A. Di- 
Nola, and J. R. Haak, J. Chem. Phys. 81, 3684 (1984) 

2® J. M. Soler, E. Artacho, J. D. Gale, A. Garcia, J. Junquera, P. Or- 
dejon, and D. Sanchez-Portal, J. Phys.: Condens. Matter 14,2745 
( 2002 ) 

20 J. P. Perdew, K. Burke, and M. Emzerhof, Phys. Rev. Lett. 77, 
3865(1996) 

2* J. Junquera, O. Paz, D. Sanchez-Portal, and E. Artacho, Phys. 
Rev. B 64, 235111 (2001) 

22 N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991) 

22 N. D. Mermin, Phys. Rev. 137, A1441 (1965) 

2^ S. C. Harvey, R. K.-Z. Tan, and T. E. Cheatham, J. Comput. 

Chem. 19, 726 (1998) 

22 J. J. Hall, Phys. Rev. 161, 756 (1967) 

20 A. Kubo, Y. Wang, C. E. Runge, T. Uchida, B. Kiefer, 
N. Nishiyama, and T. S. Duffy, Journal of Physics and Chem¬ 
istry of Solids 69, 2255 (2008) 

2' R. E. Duff and E. S. Minshall, Phys. Rev. 108, 1207 (1957) 

22 R. Hull, Properties of Crystalline Silicon, EMIS datareviews se¬ 
ries (INSPEC, the Institution of Electrical Engineers, 1999). 




